L’objectif de ce TP est d’apprendre à manipuler JDemetra+ sous R à travers le JWSACruncher et le package RJDemetra.
Pour manipuler JDemetra+ sous R il y a actuellement deux façons :
Utiliser le JWSACruncher (https://github.com/jdemetra/jwsacruncher) qui permet, à partir de la console, de mettre à jour un workspace JDemetra+ et d’exporter les résultats sans devoir ouvrir le logiciel. Pour faciliter son utilisation depuis R, le package rjwsacruncher peut être utilisé.
Utiliser le package RJDemetra qui permet d’effectuer des désaisonnalisations avec les mêmes algorithmes et paramètres que JDemetra+ et de manipuler des workspaces.
rjwsacruncherLe JWSACruncher est téléchargeable ici : https://github.com/jdemetra/jwsacruncher/releases.
Pour utiliser la dernière version il faut avoir une version de Java supérieure à la 8, si ce n’est pas le cas, il faut télécharger une version portable de Java et configurer le JWSACruncher en conséquence (https://github.com/AQLT/BCEAO_2020/wiki/Installation-et-configuration-de-JDemetra---et-du-JWSACruncher#installation-du-jwsacruncher). Ces manipulations peuvent aussi se faire à partir de rjwsacruncher :
# install.packages("rjwsacruncher") # Si pas déjà installé
library(rjwsacruncher)
# Télécharge l'archive du JWSACruncher et la met sur le D:/
download_cruncher("D:/")
# Dézipper l'archive et ensuite pour configurer avec une version portable de Java :
jwsacruncher_path <- "D:/jwsacruncher-2.2.2-bin/bin/jwsacruncher.bat" # Lien vers le fichier jwsacruncher.bat
java_path <- "D:/Java8/bin/java.exe" # Lien vers le fichier java.exe de la version portable de Java
configure_jwsacruncher(jwsacruncher_path, java_path)
Pour indiquer à rjwsacruncher où se trouve le JWSACruncher, le plus simple est de mettre à jour l’option cruncher_bin_directory :
# Chemin vers le dossier bin du JWSACruncher
options(cruncher_bin_directory =
"/Users/alainquartierlatente/Desktop/BCEAO/jwsacruncher-2.2.2-bin/bin")
getOption("cruncher_bin_directory") # Pour afficher la valeur actuelle
## [1] "/Users/alainquartierlatente/Desktop/BCEAO/jwsacruncher-2.2.2-bin/bin"
Pour lancer le JWSACruncher il faut trois fichiers :
create_param_file()) ;cruncher_bin_directory).create_param_file() pour créé un fichier de paramètres permettant de mettre à jour un workspace : - En reestimant le modèle ARIMA, les outliers et les autres paramètres du modèle de régression et en re-identifiant les outliers uniquement sur la dernière année. - En exportant la statistique M7, la statistique Q-M2 et les tests de jours ouvrables résiduels ; - En exportant La série brute, la série désaisonnalisée et la tendance (de manière verticale).
create_param_file(dir_file_param = "/Users/alainquartierlatente/Desktop/BCEAO/",
policy = "lastoutliers",
matrix_item = c("m-statistics.m7",
"m-statistics.q-m2",
"diagnostics.residual trading days tests.f-test on sa (td):2",
"diagnostics.residual trading days tests.f-test on i (td):2"),
tsmatrix_series = c("y", "sa", "t"),
csv_layout = "vtable"
)
Une fois ce fichier créé, pour lancer le JWSAcruncher il suffit d’utiliser la fonction cruncher :
cruncher(workspace = "workspace.xml",
param_file_path = "/Users/alainquartierlatente/Desktop/BCEAO/parameters.param"
)
Si vous n’avez pas de workspace vous pouvez utiliser le code suivant (qu’on expliquera plus tard) pour en générer un :
library(RJDemetra)
spec_x13 <- x13_spec(spec = "RSA5c", easter.enabled = FALSE)
sa_x13 <- x13(ipi_c_eu[, "FR"], spec = spec_x13)
spec_ts <- tramoseats_spec(spec = "RSA5")
sa_ts <- jtramoseats(ipi_c_eu[, "FR"], spec = spec_ts)
wk <- new_workspace()
new_multiprocessing(wk, "sa1")
add_sa_item(wk, "sa1", sa_x13, "X13")
add_sa_item(wk, "sa1", sa_ts, "TramoSeats")
save_workspace(wk, "workspace.xml")
Si non spécifié dans le fichier des paramètres, les résultats sont exportés dans le sous dossier “Output” du workspace (pour le workspace.xml, les résultats seront donc sous workspace/Output/). On peut aussi créer le fichier des paramètres et lancher le JWSAcruncher avec la fonction cruncher_and_param. Cette fonction permet aussi de renommer les dossiers exportées avec les noms des multi-processings utilisés dans JDemetra+ (évite d’avoir des dossiers du type SAProcessing-1) =
cruncher_and_param(
workspace = "workspace.xml",
policy = "lastoutliers",
matrix_item = c("m-statistics.m7",
"m-statistics.q-m2",
"diagnostics.residual trading days tests.f-test on sa (td):2",
"diagnostics.residual trading days tests.f-test on i (td):2"),
tsmatrix_series = c("y", "sa", "t"),
csv_layout = "vtable"
)
Pour faire de la désaisonnalisation sous R il existe plusieurs packages :
seasonal et x12 qui permettent de faire du X-13ARIMA-SEATS en utilisant les programmes du US Census Bureau
RJDemetra qui est une interface R à JDemetra+ et c’est ce package qu’on va étudier.
RJDemetra est sur le CRAN et se base sur les librairies Java de JDemetra+. Pour l’utiliser il faut avoir Java 8 ou plus. En cas de problème d’installation voir la page : https://github.com/jdemetra/rjdemetra/wiki/Installation-manual.
Le package a aussi un site web (https://jdemetra.github.io/rjdemetra/).
RJDemetra permet :
de faire des modèles RegARIMA, TRAMO-SEATS and X-13-ARIMA comme dans JDemetra+ en définissant sa propre spécification
Dans les prochains exercices, la série utilisée sera ipi_c_eu[, "FR"] qui est l’IPI français. Vous pouvez bien sûr adapter le code pour utiliser vos propres séries. Les fonctions utilisées seront x13(), x13_spec(), regarima_x13, regarima_x13_spec ou regarima.
Utiliser la spécification RSA4c pour la désaisonnalisation.
mysa, regarder les valeurs de mysa$final, mysa$final$series et mysa$final$forecasts.
mysa <- x13(ipi_c_eu[, "FR"], spec = "RSA4c")
mysa
## RegARIMA
## y = regression model + arima (2, 1, 1, 0, 1, 1)
## Log-transformation: no
## Coefficients:
## Estimate Std. Error
## Phi(1) 0.3358 0.171
## Phi(2) 0.2060 0.096
## Theta(1) -0.2450 0.173
## BTheta(1) -0.5112 0.050
##
## Estimate Std. Error
## Easter [1] -1.133 0.337
## LS (11-2008) -8.000 1.283
## LS (1-2009) -7.551 1.283
## LS (5-2008) -5.069 1.234
##
##
## Residual standard error: 1.696 on 314 degrees of freedom
## Log likelihood = -631, aic = 1280 aicc = 1281, bic(corrected for length) = 1.2
##
##
##
## Decomposition
## Monitoring and Quality Assessment Statistics:
## M stats
## M(1) 0.063
## M(2) 0.040
## M(3) 0.957
## M(4) 0.218
## M(5) 0.744
## M(6) 0.175
## M(7) 0.076
## M(8) 0.208
## M(9) 0.056
## M(10) 0.154
## M(11) 0.138
## Q 0.267
## Q-M2 0.295
##
## Final filters:
## Seasonal filter: 3x5
## Trend filter: 13 terms Henderson moving average
##
##
## Final
## Last observed values
## y sa t s i
## Jan 2017 97.4 100.5322 100.5705 -3.1322239 -0.03826496
## Feb 2017 97.5 100.2684 100.9677 -2.7684016 -0.69925325
## Mar 2017 112.0 102.6068 101.4271 9.3932126 1.17967381
## Apr 2017 103.0 101.0608 101.8755 1.9391548 -0.81462781
## May 2017 100.4 103.0181 102.2797 -2.6181258 0.73846654
## Jun 2017 111.2 102.4967 102.6808 8.7032680 -0.18410827
## Jul 2017 103.4 103.1869 103.0846 0.2130743 0.10229998
## Aug 2017 79.3 103.2404 103.5149 -23.9404354 -0.27443903
## Sep 2017 109.7 103.5634 103.9774 6.1366066 -0.41399806
## Oct 2017 114.0 106.6880 104.4467 7.3120438 2.24128542
## Nov 2017 107.7 105.5033 104.8516 2.1967210 0.65170893
## Dec 2017 101.4 104.7987 105.1825 -3.3987323 -0.38380275
##
## Forecasts:
## y_f sa_f t_f s_f i_f
## Jan 2018 102.43345 105.4677 105.3976 -3.0342512 0.070123605
## Feb 2018 102.58587 105.4460 105.5444 -2.8601225 -0.098412224
## Mar 2018 114.09123 105.7687 105.6793 8.3224915 0.089451587
## Apr 2018 108.72980 105.7481 105.8442 2.9817073 -0.096148352
## May 2018 103.49666 106.0579 106.0634 -2.5611982 -0.005554506
## Jun 2018 114.90456 106.2747 106.3311 8.6298645 -0.056376001
## Jul 2018 106.80263 106.6847 106.6316 0.1179412 0.053109111
## Aug 2018 83.03426 106.9445 106.9537 -23.9102323 -0.009215700
## Sep 2018 113.05794 106.9764 107.2581 6.0815141 -0.281665747
## Oct 2018 115.38012 108.0004 107.5237 7.3797310 0.476732719
## Nov 2018 110.13209 107.8139 107.7395 2.3182391 0.074354444
## Dec 2018 104.22637 107.6754 107.8973 -3.4490221 -0.221902905
##
##
## Diagnostics
## Relative contribution of the components to the stationary
## portion of the variance in the original series,
## after the removal of the long term trend
## Trend computed by Hodrick-Prescott filter (cycle length = 8.0 years)
## Component
## Cycle 1.652
## Seasonal 49.432
## Irregular 0.441
## TD & Hol. 0.022
## Others 49.599
## Total 101.145
##
## Combined test in the entire series
## Non parametric tests for stable seasonality
## P.value
## Kruskall-Wallis test 0.000
## Test for the presence of seasonality assuming stability 0.000
## Evolutive seasonality test 0.015
##
## Identifiable seasonality present
##
## Residual seasonality tests
## P.value
## qs test on sa 1.000
## qs test on i 1.000
## f-test on sa (seasonal dummies) 0.998
## f-test on i (seasonal dummies) 0.983
## Residual seasonality (entire series) 0.995
## Residual seasonality (last 3 years) 0.936
## f-test on sa (td) 0.001
## f-test on i (td) 0.003
##
##
## Additional output variables
y <- mysa$final$series[,"y"]
# De façon équivalente :
y <- get_ts(mysa)
sa <- mysa$final$series[,"sa"]
plot(y)
lines(sa, col = "red")
# ou on peut directement utiliser les fonctions de RJDemetra :
plot(mysa, first_date = 2000, #Pour n'afficher le graphique qu'à partir de 200
type_chart = "sa-trend" # Pour faire le graphique avec y, sa et tendance
)
spec_sans_easter <- x13_spec(mysa,
easter.enabled = FALSE)
mysa2 <- x13(ipi_c_eu[, "FR"], spec_sans_easter)
mysa2$regarima
## y = regression model + arima (0, 1, 1, 0, 1, 1)
## Log-transformation: no
## Coefficients:
## Estimate Std. Error
## Theta(1) -0.5587 0.047
## BTheta(1) -0.5279 0.049
##
## Estimate Std. Error
## LS (11-2008) -8.496 1.340
## LS (1-2009) -7.425 1.340
## LS (5-2008) -5.390 1.274
##
##
## Residual standard error: 1.742 on 317 degrees of freedom
## Log likelihood = -639.7, aic = 1291 aicc = 1292, bic(corrected for length) = 1.2
summary().
summary(mysa2$regarima)
## y = regression model + arima (0, 1, 1, 0, 1, 1)
##
## Model: RegARIMA - X13
## Estimation span: from 1-1990 to 12-2017
## Log-transformation: no
## Regression model: no mean, no trading days effect, no leap year effect, no Easter effect, outliers(3)
##
## Coefficients:
## ARIMA:
## Estimate Std. Error T-stat Pr(>|t|)
## Theta(1) -0.55869 0.04657 -12.00 <2e-16 ***
## BTheta(1) -0.52795 0.04867 -10.85 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Regression model:
## Estimate Std. Error T-stat Pr(>|t|)
## LS (11-2008) -8.496 1.340 -6.339 7.77e-10 ***
## LS (1-2009) -7.425 1.340 -5.541 6.25e-08 ***
## LS (5-2008) -5.390 1.274 -4.232 3.02e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 1.742 on 317 degrees of freedom
## Log likelihood = -639.7, aic = 1291, aicc = 1292, bic(corrected for length) = 1.2
Dans cette partie nous allons créer un workspace depuis R. Pour cela les fonctions qui peuvent être utilisées sont new_workspace(), load_workspace(), new_multiprocessing(), add_sa_item(), save_workspace(), compute(), get_object(), get_name(), get_ts() ou count.
wk <- new_workspace()
new_multiprocessing(wk, "MP-1")
add_sa_item(wk, "MP-1", mysa, "X13 avec Pâques")
add_sa_item(wk, "MP-1", mysa2, "X13 sans Pâques")
add_sa_item(wk, "MP-1", tramoseats(ipi_c_eu[, "FR"]), "TRAMO-SEATS")
save_workspace(wk, "mon_premier_workspace.xml")
wk <- load_workspace("mon_premier_workspace.xml")
compute(wk)
count(wk) # Nombre de multiprocessing
## [1] 1
mp1 <- get_object(wk, 1) # Le premier multiprocessing
get_name(mp1)
## [1] "MP-1"
count(mp1)
## [1] 3
all_y <- get_ts(mp1) # toutes les séries brutes
model2 <- get_object(mp1, 2) # On récupère l'objet associé au 2ème modèle
get_model(model2, wk)
## RegARIMA
## y = regression model + arima (0, 1, 1, 0, 1, 1)
## Log-transformation: no
## Coefficients:
## Estimate Std. Error
## Theta(1) -0.5587 0.047
## BTheta(1) -0.5279 0.049
##
## Estimate Std. Error
## LS (11-2008) -8.496 1.340
## LS (1-2009) -7.425 1.340
## LS (5-2008) -5.390 1.274
##
##
## Residual standard error: 1.742 on 317 degrees of freedom
## Log likelihood = -639.7, aic = 1291 aicc = 1292, bic(corrected for length) = 1.2
##
##
##
## Decomposition
## Monitoring and Quality Assessment Statistics:
## M stats
## M(1) 0.074
## M(2) 0.053
## M(3) 1.032
## M(4) 0.520
## M(5) 0.766
## M(6) 0.200
## M(7) 0.077
## M(8) 0.218
## M(9) 0.056
## M(10) 0.170
## M(11) 0.153
## Q 0.308
## Q-M2 0.340
##
## Final filters:
## Seasonal filter: 3x5
## Trend filter: 13-Henderson
##
##
## Final
## Last observed values
## y sa t s i
## Jan 2017 97.4 100.5913 100.6459 -3.19125580 -0.05468575
## Feb 2017 97.5 100.3932 101.0670 -2.89322193 -0.67373301
## Mar 2017 112.0 102.6491 101.5555 9.35088124 1.09357231
## Apr 2017 103.0 100.6322 102.0059 2.36777259 -1.37370330
## May 2017 100.4 103.1260 102.3959 -2.72597512 0.73006300
## Jun 2017 111.2 102.5677 102.7665 8.63229838 -0.19884777
## Jul 2017 103.4 103.3222 103.1377 0.07777798 0.18452040
## Aug 2017 79.3 103.0579 103.5427 -23.75790898 -0.48477871
## Sep 2017 109.7 103.6142 103.9860 6.08584154 -0.37187192
## Oct 2017 114.0 106.7526 104.4289 7.24736328 2.32369073
## Nov 2017 107.7 105.4763 104.7813 2.22365154 0.69502454
## Dec 2017 101.4 104.7798 105.0285 -3.37983030 -0.24869062
##
## Forecasts:
## y_f sa_f t_f s_f i_f
## Jan 2018 101.94580 105.0435 105.1462 -3.09766844 -0.10273186
## Feb 2018 102.23454 105.2355 105.2166 -3.00092915 0.01887267
## Mar 2018 114.62722 105.2221 105.3128 9.40512478 -0.09068921
## Apr 2018 107.62003 105.2912 105.4822 2.32879858 -0.19099122
## May 2018 103.14889 105.8337 105.7299 -2.68484181 0.10379097
## Jun 2018 114.60285 106.0633 106.0258 8.53958658 0.03741542
## Jul 2018 106.49502 106.5160 106.3350 -0.02101655 0.18103143
## Aug 2018 82.69293 106.3992 106.6426 -23.70626502 -0.24336723
## Sep 2018 112.73049 106.7061 106.9199 6.02439223 -0.21377196
## Oct 2018 114.98821 107.6746 107.1619 7.31365925 0.51266621
## Nov 2018 109.74540 107.3822 107.3610 2.36323681 0.02119759
## Dec 2018 103.87554 107.3301 107.5102 -3.45455021 -0.18015246
##
##
## Diagnostics
## Relative contribution of the components to the stationary
## portion of the variance in the original series,
## after the removal of the long term trend
## Trend computed by Hodrick-Prescott filter (cycle length = 8.0 years)
## Component
## Cycle 1.622
## Seasonal 48.384
## Irregular 0.439
## TD & Hol. 0.000
## Others 51.257
## Total 101.703
##
## Combined test in the entire series
## Non parametric tests for stable seasonality
## P.value
## Kruskall-Wallis test 0.000
## Test for the presence of seasonality assuming stability 0.000
## Evolutive seasonality test 0.012
##
## Identifiable seasonality present
##
## Residual seasonality tests
## P.value
## qs test on sa 1.000
## qs test on i 1.000
## f-test on sa (seasonal dummies) 0.982
## f-test on i (seasonal dummies) 0.913
## Residual seasonality (entire series) 0.987
## Residual seasonality (last 3 years) 0.965
## f-test on sa (td) 0.001
## f-test on i (td) 0.004
##
##
## Additional output variables
L’objectif de cette partie est de manipuler la fonction jx13() pour gagner en temps de calcul.
jx13() et la spécification sans effet graduel de pâques calculée dans une des sections précédentes.
myjsa <- jx13(ipi_c_eu[, "FR"], spec_sans_easter)
get_indicators(myjsa, "sa")
## $sa
## Jan Feb Mar Apr May Jun Jul
## 1990 93.38868 95.34727 94.66473 94.22687 94.82241 93.53732 93.40578
## 1991 93.72301 92.42676 92.80794 92.30200 91.03344 93.23475 91.13076
## 1992 92.15009 91.99782 92.60121 91.97359 91.96023 91.47582 90.50553
## 1993 88.02599 87.55694 86.45328 86.81808 86.43424 85.47874 86.91383
## 1994 87.75673 87.50119 87.37781 89.19783 89.85261 90.34048 91.13467
## 1995 93.46692 94.21319 93.60721 93.77304 92.73745 93.41740 93.79053
## 1996 93.59156 92.70057 93.67606 93.26070 94.30829 93.77376 92.05111
## 1997 93.96521 95.72488 96.62393 99.37782 96.94193 97.79877 97.75854
## 1998 102.85584 103.23915 101.47488 102.82653 103.78015 103.40683 104.40880
## 1999 104.25168 102.83976 103.28376 104.10695 104.43591 105.69456 106.11420
## 2000 108.58303 109.00558 109.65561 109.74390 111.46942 110.65188 111.41231
## 2001 112.60684 113.09267 113.85094 110.95996 112.22599 112.84385 110.78450
## 2002 110.35644 110.69634 111.12976 111.49529 110.61059 110.38909 109.52486
## 2003 109.12552 109.37528 109.78055 109.50054 107.06669 106.38891 108.69508
## 2004 109.44578 110.82064 109.82648 110.44698 109.84808 111.54939 111.86543
## 2005 112.77430 110.65979 108.13258 112.06755 110.30359 110.01187 109.93903
## 2006 110.89634 110.26646 112.60474 111.41607 113.71308 114.33799 111.84987
## 2007 113.05741 113.91173 114.91586 113.39637 114.95875 114.75922 114.18828
## 2008 114.80640 115.84724 114.35443 116.07897 110.46594 110.43233 111.15412
## 2009 95.00383 93.74985 91.55923 92.31293 93.44255 94.57255 94.87634
## 2010 96.87329 96.54526 98.33378 98.60002 99.06908 98.76055 99.51520
## 2011 103.65514 104.89709 103.47376 101.89841 105.21026 101.74471 102.74365
## 2012 101.21018 100.30565 101.57684 100.08670 99.46907 99.57143 100.47658
## 2013 97.82924 98.48466 97.58849 99.93840 99.46724 99.16283 99.04460
## 2014 98.54269 99.74664 98.59520 99.00308 96.36291 98.73647 99.48021
## 2015 98.77168 99.23059 99.83210 99.22749 99.28563 101.15290 99.23240
## 2016 101.84644 100.65843 98.21419 100.89193 99.64342 99.03778 99.32238
## 2017 100.59126 100.39322 102.64912 100.63223 103.12598 102.56770 103.32222
## Aug Sep Oct Nov Dec
## 1990 92.39240 93.29289 93.49084 91.67249 91.82757
## 1991 91.96046 92.54231 91.53246 93.18318 91.79487
## 1992 90.92790 89.26406 89.47289 89.61672 88.84348
## 1993 86.93046 86.49104 86.85764 85.55160 86.30094
## 1994 91.46393 91.12375 92.09164 94.18695 93.99897
## 1995 92.84402 93.94477 92.83151 92.24910 95.61548
## 1996 95.11107 94.28116 93.11313 94.58966 94.76052
## 1997 101.93484 99.79300 102.42193 100.03530 103.22334
## 1998 102.38845 103.50021 102.99348 103.64748 103.45541
## 1999 104.41596 107.16166 108.37309 108.73583 108.91223
## 2000 110.23237 111.03054 111.72169 112.32448 115.44455
## 2001 114.03205 110.97730 110.97886 109.82051 109.97363
## 2002 112.17047 109.81120 108.54410 111.44029 107.78987
## 2003 108.68313 108.03508 109.79495 109.41776 109.70771
## 2004 108.38866 111.75950 111.89501 110.76189 110.45906
## 2005 109.88989 113.00823 108.49387 114.09147 112.07881
## 2006 111.99666 113.86033 112.29181 112.58976 114.25597
## 2007 114.58119 113.05487 115.63052 114.03865 113.92006
## 2008 109.95485 109.56795 107.66524 101.11017 101.09179
## 2009 98.83345 97.12055 96.69080 97.98809 96.41463
## 2010 99.65046 100.12052 100.51787 101.08592 101.50428
## 2011 102.42480 101.51561 102.46471 104.61316 101.87293
## 2012 102.74682 99.97581 98.14837 97.93569 99.04334
## 2013 99.01646 98.11041 99.24240 99.64230 98.36856
## 2014 98.59309 99.15174 98.43846 96.92921 99.96872
## 2015 100.91215 100.83953 101.07252 99.98627 100.50965
## 2016 101.35740 100.06667 99.36811 102.14185 100.42113
## 2017 103.05791 103.61416 106.75264 105.47635 104.77983
window(time(ipi_c_eu[, "FR"]), start = 2005).
dates <- window(time(ipi_c_eu[, "FR"]), start = 2005)
estimations <- sapply(dates, function(last_date_estimation){
myjsa <- jx13(window(ipi_c_eu[, "FR"], end = last_date_estimation), spec_sans_easter)
sa <- get_indicators(myjsa, "sa")$sa
window(sa, start = 2005, end = 2005) # Pour ne récupérer que la valeur en 2005
})
estimations <- ts(estimations, start = 2005, frequency = 12)
plot(estimations)